J. Chem. Phys., submitted 11/12/03. 



An Empirical Charge Transfer Potential with 
Correct Dissociation Limits 

Steven M. Valone 

Materials Science and Technology Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545 and 
Department of Physics and Astronomy, University of New Mexico, Albuquerque, New Mexico 87131 

Susan R. Atlas 

Center for Advanced Studies and Department of Physics and Astronomy, 
University of New Mexico, Albuquerque, New Mexico 87131 

(Dated: February 2, 2008) 

The empirical valence bond (EVB) method [J. Chem. Phys. 52, 1262 (1970)] has always embodied 
charge transfer processes. The mechanism of that behavior is examined here and recast for use as 
a new empirical potential energy surface for large-scale simulations. A two-state model is explored. 
The main features of the model are: (1) Explicit decomposition of the total system electron density is 
invoked; (2) The charge is defined through the density decomposition into constituent contributions; 
(3) The charge transfer behavior is controlled through the resonance energy matrix elements which 
cannot be ignored; and (4) A reference-state approach, similar in spirit to the EVB method, is 
used to define the resonance state energy contributions in terms of "knowable" quantities. With 
equal validity, the new potential energy can be expressed as a nonthermal ensemble average with 
a nonlinear but analytical charge dependence in the occupation number. Dissociation to neutral 
species for a gas-phase process is preserved. A variant of constrained search density functional 
theory is advocated as the preferred way to define an energy for a given charge. 
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I. INTRODUCTION 

Charge transfer is ubiquitous in physical processes af- 
fecting biological, chemical, and materials systems. The 
representation of charge transfer is of intense current in- 
terest throughout the physical sciences. A powerful con- 
cept in both modeling and understanding how charges re- 
distribute themselves during a physical process is chemi- 
cal potential equalization. To apply chemical potential 
equalization successfully, it is essential to use a charge- 
dependent energy model which behaves correctly for all 
configurations encountered in the process of interest. 

Consider a diatomic molecule AB. Atom A is assumed 
to be more electropositive than atom B. We are interested 
in charge disproportionation reactions typified by 

A°B" ^ A+^B-' (1) 

and 

A2+B2 ^2A+«B-« , (2) 

where the charge q is not necessarily an integer. Charge 
disproportionation and transfer reactions are impor- 
tant in electron transfer in biophysical systems, en- 
zyme catalysis reactions, ^"'^^ and such processes as 
electron-hole production and recombination in organic 
semiconductors.^^ For dimers (B = A), Eq. (1) cor- 
responds to broken charge-symmetry states. Further- 
more, Eqs. (1) and (2) have been the prototype reac- 
tions for a wide variety of theoretical studies of chemical 
bonding.^'^^"^° 



The most prevalent model potential for charged species 

is the quadratic expansion of the classical electrostatic 
potential,'''6'2i"23 

E « EoiR) + Ei{R) q + 1/2 {E2{R) - V{R)) . (3) 

The Ei{R) are expansion coefficients which are atom- 
type dependent and V{R) is an interionic potential repre- 
senting the classical electrostatic contribution. In general 
all of these functions depend on the separation R. The 
charge q is usually, but not always, independent of R. 
V{R) approaches a \/R dependence at sufficiently large 
R. Eo{R) contains the charge- independent short-range 
and dispersion interactions between the atoms. The 
other expansion coefficients are frequently interpreted in 
terms of physical quantities such as chemical potential 
and hardness. ^^^'^'^"'^"'^ The quadratic form has the virtue 
of simplicity and works adequately when the range of R 
is small enough to prevent q from changing appreciably. 

If Eq. (3) is used at all separations, the atoms in the 
AB molecule will remain ionic even at separations where 
they are supposed to return to neutral states. For those 
situations where charge transfer does occur, alternative 
functional forms of charge dependence need to be in- 
voked. For instance, Alavi et al.^^ and Grochowski and 
coworkers^^ use phenomenological switching functions to 
effect charge transfer. At very large separations, though, 
it is known that the charge dependence becomes piece- 
wise linear (Fig. 1).2>19^33 

Morales and Martinez'^^ (hereafter referred to as MM) 
conclude that the "grand canonical" (GC) approach^ 
cannot describe a realistic charge transfer process. In the 
GC approach each atom is considered to be an open sys- 
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FIG. 1: Comparison of piecewise linear and quadratic mod- 
els of the charge dependence of the energy for a diatomic. 
The solid line corresponds to the energy of the isolated atoms 
as a function of acquired charge q (measured with respect 
to charge transfer from the cation). The dashed curve cor- 
responds to a quadratic representation of the charge depen- 
dence. The energy minimizes at g = for the piecewise linear 
model, but not for the quadratic model. At large separa- 
tions, the piecewise linear model is correct. Traditionally the 
quadratic model is considered to be more correct at finite 
separations. 



tcm with respect to exchange of energy and numbers of 
electrons. They use a 3-state model with integer charge 
resonances representing the states of the ensemble for 
each atom. For atom A for instance, the GC energy is 
expressed in the form of an ensemble average 



E' 



GC 



(4) 



where are the occupation numbers and Ea" are the 
energies for the integer charge species a — -f , 0, and — . 
The LHA" depend on q and are equal to or greater than 
zero. (See Refs. 26, 34, and 35 for the detailed expres- 
sions.) MM show that the GC energy as expressed by 
Eq. (4) minimizes to integer charges, and never to frac- 
tional charges. When the two atoms are well separated 
and only weakly interacting, the covalent and ionic reso- 
nance states can be regarded as being close to eigenfunc- 
tions and the mixing terms between the states can be 
ignored. The nonlinearities in the occupation numbers 
as functions of charge introduced through temperature 
at physically reasonable temperatures are not sufficient 
to produce states with fractional charges. Thus charge 
transfer as represented by an ensemble of weakly inter- 
acting integer-charge states is too simplistic to describe 
fractionally charged states. 

MM also examine a 3-state valence bond (VB) ap- 
proach usqing the same resonance states as in the GC 



approach. This is equivalent to a K ensemble where 
the states are allowed to interact. The VB approach 
is able to represent charge transfer processes with frac- 
tional charges. However, the representation of the en- 
ergy depends on resonance matrix elements. The issue 
becomes one of defining sets of coefficients for the state 
wavcfunctions that allow one to recover quadratic and 
GC-like energy expressions. Using Coulson charges^^"'^* 
as an approximation and a maximum entropy valence 
bond approach (MEVB), MM derive energy expressions 
as functions of q with the forms of Eqs. (3) and (4). In 
the MEVB approach, decomposition into atomic contri- 
butions is accomplished through an examination of the 
various matrix elements appearing in the total energy 
expression. However, MM were unable to find a gen- 
eral expression linking the GC and classical electrostatic 
forms. 

The three main difficulties in deriving a charge depen- 
dent potential are the same ones facing MM and others. 
First is to define the charge, second is to evaluate or elim- 
inate the resonance energy matrix elements, and third is 
to define the energy for a particular value of q. First, 
to define g, we invoke a density- functional-theory'^^"'*^ 
(DFT) motivated atom-in-moleculc"'^*'^^'''^ pseudo-atom 
concept within the context of the EVB approach. We 
assume the availability of a practical density decomposi- 
tion strategy^*'^'^"** to define the pseudo-atom densities. 
The charge is defined as an average over the difference 
between two pseudo-atom densities. No restriction to 
Coulson'^^ or other definitions of charge^^'^" is necessary. 
The dependence of the EVB wavefunction on q is then de- 
duced, which in turn yields an energy for arbitrary not 
just the optimum q. Second, to evaluate resonance en- 
ergy matrix elements, we retain the VB approach of MM, 
but use an empirical valence bond (EVB)-^^'^*'^^ strategy, 
rather than explicitly evaluating a model Hamiltonian or 
using the MEVB averaging procedure of MM. Specifi- 
cally, a reference energy is separated out for a fixed value 
of q. To define the energy for an atom in the molecule, we 
require consistency between the EVB wavefunction en- 
ergy and DFT energy for each atom. In so doing we are 
able to cast the energy of each atom into a form suitable 
for constructing empirical potential energy surfaces that 
could be used in simulations of larger systems. Third, 
the typical situation is that there are many wavefuctions 
and many electron densities that are consistent with a 
particular q. To define a unique energy from among the 
possible choices, we adhere to constrained search density 
functional theory (CS-DFT)''*'^^ rather than appealing 
to a maximum entropy principle. 

To the best of our knowledge, the EVB method has 
not previously been combined with an atom-in-molecule 
approach. von Szentpaly and coworkers define atomic 
charges, but their definition is implicitly limited to a 
quadratic dependence. ^^"^^ Grochowski and coworkers^^ 
use EVB parameterized by ah initio calculations for key 
molecular fragments and for defining the charges. The 
charges are defined with a phenomenological spatial de- 
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pendence that docs maintain correct dissociation to neu- 
trals. Furthermore, the charge dependence of the final 
potential is purely electrostatic. 

The inconsistency between the GC and classical elec- 
trostatic forms is explained in Cioslowski and Stefanov^^ 
using a different definition of q, which is based on the 
total system wavefunction. The charge is expressed as 
a perturbation on the molecular Hamiltonian H. Com- 
puted atom-in-molecule equilibrium charges are used in 
subsequent "charge-constrained" calculations to study 
the energy and electronegativity dependences about the 
ground-state atom-in-molecule charges. Nalewajski^"^ 
provides a simplified rendition of Cioslowski and Ste- 
fanov, but only for orthogonal resonance states. While 
physically correct, these efforts have not been cast in the 
form of general purpose potential energy surfaces. In par- 
ticular, neither attempts to resolve the resonance energy 
issue in a way that is tractable for large-scale simulations. 

Here wc derive a charge-dependent empirical poten- 
tial which faithfully represents charge transfer as a func- 
tion of separation between atom A and an entity B. In 
the simplest case, entity B is another atom. More gen- 
erally, B represents a collective embedding environment 
or reservior.^'^^ The charge may be fractional. Our pri- 
mary interest is in deriving a general functional form with 
correct physical and chemical behavior at all interaction 
strengths rather than providing an exact treatment of 
particular terms or systems. The three difficulties just 
outlined are addressed. We then analyze these potentials 
in special limiting cases and in the light of the results of 
Cioslowski and Stefanov^^ and Nalewajski.'^^ Finally, we 
derive general models for both pair potentials and atom- 
in-molecule energies. 



II. THE EMPIRICAL VALENCE BOND 
REPRESENTATION 



The EVB method is a much more general technique 
than described here. Reviews of EVB are available from 
Warshel et. al.^^ Here we confine our discussion to a 2- 
state model of the AB molecule. There are fixed covalcnt 
and ionic resonance states represented by wavefunctions 
tpc and ipi, respectively. ■^■^'^■^'^'''^^ No assumption is made 
about the quality of these wavefunctions. However, ipc 
retains A'ao electrons on atom A and N^o electrons on 
atom B. Atoms A and B are neutral in the covalent res- 
onance state. Consequently N/^o and A^b^ are equal to 
their respective nuclear charges, and Z^. Similarly, 
tpi retains N^o — 1 electrons on atom A and Nqo + 1 
electrons on atom B. In each resonance state, the total 
number of electrons is A. tpc and ipi are assumed to be 
normalized to unity. The wavefunction of the system tp 
is the combination^^"^^'^^"^^ 



where 7 determines the relative ionic character. By nor- 
malization. 



where 



l/c2 = 1 + 275^+7' , 



Sci = (V'clV'i) 



(6) 



(7) 



For the AB molecule with Hamiltonian H, the mixed- 
state energy takes the form 

E = c2(i/,, + 27H„+7'i?n) , (8) 



1 + 275',, + 72 



(9) 



where Hcc, Ha and Hd are the associated energy matrix 
elements. 



and 



(10) 

(11) 

(12) 



Minimizing E with respect to 7 gives the optimized 
values, 



7opt — 



1 ± y/1 + tec Ci 



(13) 



where Ecc = 2 (iJcj - SciHcc)/{Hu ~ Hcc) and en = 
2 {Hc^-Sc^ Hu)/iHu-Hcc). Thc ± signs in Eq. (13) cor- 
respond to ground (gs) and excited (xs) states, whose en- 
ergies are designated as Egg and E^g, respectively. From 
Eq. (13), one can see that the off-diagonal matrix el- 
ements, Hci and Sd, control charge transfer in EVB. 
When Scr = 0, ecc = = e = 2H„/{Hu - Hcc)- De- 
pending on the root, 7opt (— root) or l/7opt (+ root) 
varies from —1 to +1. Note that as Ha — Hcc goes 
to zero, that is, as the covalent and ionic curves cross, 
charge transfer or the Coulson-Fischer transition^^ be- 
comes more abrupt. If ipc and ipi are nondegenerate 
eigenfunctions of H , then Hd = Sd = 0. There is either 
no charge transfer or there is complete charge transfer. 
In either case, the states do not mix. 

The coefficient 7 governs the ionic contribution to ip. 
As 7 increases, ip migrates toward more complete charge 
transfer. Even neglecting overlap, ^^'^^ the ionic strength 
7 interpolates between covalent and ionic states in a phys- 
ically reasonable way.^^ The state-to-statc interpolating 
behavior of 7 in Eq. (13) is the essential behavior that we 
wish to emulate in developing a more broadly applicable 
charge dependent potential. 

Later in the paper, two other relationships will be- 
come useful, which we provide here. First, in the spirit 
of EVB, ^2 one always wants to know the resonance en- 
ergy in the terms of the ground-state i?gs, Sd^ Hcc, and 
Hji. That is, 



■0 = c(-0c + 



(5) 



Hr 



Er Sc^ ± ^{Hcc-Er)(Hu-Er) , (14) 
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where r is either "gs" or "xs". The second relationship The total electronic energy expression analogous to 

expresses 7opt as a function of these same parameters. Eq. (17) is 
This is done by substituting Eq. (14) into the expressions 

for and en. The general expression is E = a^c Ucc + 2 ad Hd + OLa Ha . (18) 



{Er — Haa) Sd ± \/ {Hcc — ET){Hii — Er) 



(15) 

where aa is either "cc" or "n" and either the ground or 
excited state is selected. 



III. CHARGE IN THE 2-STATE EVB MODEL 

We begin with the question of how to define charge in 
our model. There is no unique definition of the charge on 
an atom in a molecule. To assign charges to individual 
atoms, we assume that it is possible to decompose the 
total electron density into pseudo-atom densities. How- 
ever, we do not need to specify a particular decomposi- 
tion procedure at this time. We only need to know that 
some procedure is available. 



Eq. (18) corresponds to a non-diagonal representation of 
the energy. Its advantage is that the component contri- 
butions of the essential states that are thought to rep- 
resent the physical system are delineated. A diagonal 
representation of T when state mixing is important can 
be devised by diagonalizing H. For the purposes of 
extending the present treatment to finite temperature 
ensembles, ^'^^'"^^'^^ this particular diagonalization of V 
would be useful. 

The A'^-electron densities of interest are p, pcc, Pci, and 
Pa. They bear the usual relationships to the respective 
A^-electron density matrices, V„t'- 



(19) 



where rjv_i are the 3(iV-l) dimensional spatial electronic 
coordinates and ctt is either "cc" , "ci" , "ii" , or no sub- 
script, p. Pec, and pa are normalized to N. The relation- 
ship of central interest is 



A. Definition of Charge 

To define the charge, it is convenient to use iV-electron 
density matrices and 1-electron densities. In terms of 
density matrix language, the state of AB is represented 
as 

r(r^,rjv) = V(rW)^(rjv) , (16) 

where r^r are the 3N dimensional spatial electronic coor- 
dinates for the full AB system and ?/' is given in Eq. (5). 
Spin is ignored at this point and, for simplicity, the ma- 
trix elements are assumed to be real. Eq. (16) corre- 
sponds to the pure-state representation of the density 
matrix for the 2-state model. We can expand Eq. (16) 
in terms of the covalent and ionic states, resulting in the 
relationship 

r = ttcc Tec + 2 Uci Td + a^i Tii , (17) 

where 



r 

^ cc 






r . 

Cl 


= i>cipt , 




Tu 


= , 




Otcc 


= 1/(1 + 275,,- 




ad 


= 7/(1 + 275,,- 


f7^) 



and 

au=lVil + 2-fSc^+J^) . 



P = acc Pec + 2 act Pa + au Pu 

_ Pcc + 'il Pci + I^PU . 

" l + 2JSc^+7' ■ ^ ' 

Recall the assumption that the covalent and ionic state 
wavefunctions are given and fixed. Thus the total density 
p is determined solely by the value of 7. 

The energies of pcc and pu are well-defined in a con- 
ventional DFT sense. ^^'"'^ However, the "interference 
density" Pd does not have a well-defined energy in 
DFT. Nevertheless, its energy may be inferred from the 
energies of p, pcc, and pu, as we will show below. 

Next we assume that all of the p^r can be decom- 
posed into corresponding pseudo-atom densities, p*^ ^ 
and p*^B.^'^ The pseudo-atom densities p\ and pg inte- 
grate to non-integer values, and Nq, whereas P*cA' 
Pii A' Pcc Bi ^^'^ Pii B Constrained to integrate to in- 
teger numbers of electrons. We use asterisks throughout 
to indicate atom-in-molecule quantities. 

With these definitions in place, we define q from either 
pair of total and covalent pseudo-atom densities, which 
others have sometimes referred to as pseudo-atom distor- 
tion densities. ^^'^"^ We choose atom A: 

q = J dr(p:,^A(r)-pI(r)) (21) 
= Na-NX. (22) 

The density decompositions must be constrained to 
yield the correct number of electrons prescribed by 
Eq. (22).^^"^^ Finally, we note that the present defini- 
tion of q falls into Truhlar's Class II category. ^^'^^ 
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Now we want to eliminate 7 in favor of q. Assuming 
that a component definition of p\{r) based on Eq. (20) 
is possible, 



7 



2jSN* 



1 + 27^^ + 72 



where 



(r)-(7V^o/iV)p,,(r)) 



(23) 



(24) 



A^^o is the number of electrons on atom A when it is in 
a neutral state, and the relationship N Sd = J drpcj(r) 
has been used. The quantity SN*^ ^ is determined by the 
difference between p*^ ^ and the atom A component of 
the decomposition of pd with locally unbiased, electron- 
number decomposition.®'^ Clearly, different density de- 
composition strategies will yield somewhat different val- 
ues of SN*- ^. One possibility for determining SN*- ^ is 
to require consistency with the ground-state value of 7. 
For instance, if g = is optimum as for a dimer, 7 equals 
either or 2SN*^j^. For Eq. (23) to be applied suc- 
cessfully, 2 SN*^ ^ would have to correspond to a lower 
energy state than 7 = and the value of 7 would have to 
be determined from a separate calculation, such as rep- 
resented by Eq. (13). In such an approach, one would 
be effectively modeling 6N*^ ^ via a correspondence with 
the Coulson-Fischer transition. Alternatively, by anal- 
ogy with yO, we assume that single-particle determinants 
(e.g. Kohn-Sham determinants^^) can be calculated for 
the ^."^^'^^^^^ In a Kohn-Sham based approach, one 
would be effectively estimating p*^ ^ from the overlaps 
of these determinants. Additional ambiguity in 5N*^ ^ 
arises from the choices for ipc and ipi. These ambiguities 
lie behind the designation of the present approach as an 
empirical one. However, these ambiguities can be miti- 
gated by using a reference state as discussed in the next 
Subsection. To that end, it is useful to invert Eq. (23) 
so that 5N*^ ^ becomes a function of q and 7. That re- 
lationship is 



SN, 



ci,A 



(l-g)7'-25,,7-g 
27 



(25) 



It should be understood that, in the limit that 7^0, 
'5^c* A ^ also. 



B. Constructing Pair Potentials 



charge on separation for some reference state. Some em- 
pirical potentials such as EVB and the modified embed- 
ded atom method (MEAM)®^"^-'^ utilize reference states 
as a model calibration method. The methods of Mc- 
Donald and coworkers, 2^ McCammon, Grochowski, and 
coworkers, and Broughton and Mehl^^ effectively make 
q bond-length dependent. Eq. (23) provides a basis in 
EVB theory for their phenomenological charge transfer 
switching functions. 

A more attractive option is to solve for 7 — 7(g), 



7 



((57V* , A + q Sc^) ± ^ /(5iV* ,A + qS.,)^ + q{l- q) 



1-g 



(26) 

The coefficient 7 determines the strength of the contri- 
bution of ipi to -0- Eq- (26) states how the charge gov- 
erns that strength. This expression is consistent with the 
results of Cioslowski and Stefanov,^^ which are derived 
from a perturbative technique. As noted previously, even 
if q = 0, 7 equals either or 25N*^^. This is because 
the EVB model describes state mixing even when there is 
no charge transfer. For example, in the II2 molecule, the 
covalent and ionic wavefunctions mix at all finite separa- 
tions, but the ground state never involves charge transfer. 
Significantly, this formula also describes deviations from 
the ground-state charge. 

Eq. (26) can be substituted into Eq. (9), and the vari- 
ational procedure repeated. The result is the same as 
solving for q in terms of the resonance and overlap ma- 
trix elements obtained by equating Eq. (26) and (13). 

The EVB strategy is to use experimental information 
to eliminate the resonance energy.^^'^^ Here the analo- 
gous procedure is to choose a particular value oi q = qa 
and solve for Hd in terms of E(qo) and the diagonal 
matrix elements for each R. (The R dependence is sup- 
pressed.) The result is 



2aci{qo) 



(27) 



To construct a potential energy surface for AB, one 
option is to use Eq. (23) to model the dependence of the 



Substituting Eq. (27) into Eq. (18), the total energy for 
arbitrary q > has the form 



J 



E{q) = {aci{q)/aci{qo)) E{qo) + {acc{q) - acc{qo) ctcAq) / aci{qo j) Hcc + (auiq) - au{qo) aci{q) / aci{qo)) Hu ■ (28) 
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This form might be used in heu of classical electro- 
static potentials like Eq. (3) that have been in common 
use. It has the structure of an ensemble average,^ but 
the coefficients of H^c and Ha are not necessarily posi- 
tive semidefinite. Eq. (28) is constructed to possess the 
proper changes in atomic charges in the limit of molecu- 
lar dissociation. The construction of a proper ensemble 
representation is discussed in Section 5. A. 

The foundations for the quadratic dependence of the 
energy on charge, such as Eq. (3), must stem from 
Eqs. (26) and (9). However, even this simplest exam- 
ple of charge transfer has a considerably more complex 
dependence than quadratic. That dependence is clearly 
carried through the overlap contributions. 

Once go has been chosen, the procedure for determin- 
ing a point on the potential energy surface for arbitrary 
values of R and q is as follows. From some other source(s) 
of information, one must have available five reference val- 
ues: E{qQ), its associated ionicity 70, i?cc. Ha, and Sd- 
For a chosen R, one first evaluates 

(a) Qfcc, etc. from Eq. (17); 

(b) Hci from Eq. (27); and 

(c) SN:,^j, from Eq. (25). 

Then, for each q of interest, one evaluates 

(d) 7(g) from Eq. (26) using values from Steps (a), (b), 
and (c); and 

(e) E{q) from Eq. (8). 

The procedure is repeated for each value of R of interest. 

An important question is what to choose for go- One 
possible choice for go is the optimum value gopt. How- 
ever, gopt depends on Hd through Eqs. (13) and (23). 
This variant is equivalent to using the ground-state wave- 
function as one of the basis functions in the original for- 
mulation of the problem. Eq. (28) then characterizes 
deviations of the energy from the ground-state energy, 
^'('Zopt), as a function of g. Note that by choosing E{qo) 
to be consistent with E{qopt)i this variant and Eq. (28) 
will be identical. 

When go = gopt , the procedure for determining a point 
on the potential energy surface is substantially the same 
as the first. Again five reference values are needed, except 
that knowledge of gopt replaces knowledge of 70. Rela- 
tive to the first procedure, Steps (a) and (b) become to 
evaluate 

(a) Hci from Eq. (14); and 

(b) 7opt(gopt) from Eq. (13). 

All of the other steps remain the same. 

One can extend this model to g between —1 and 0. The 
entire procedure with atom A assumed to be anionic in 
is repeated. In Eq. (22), for g < 0, one simply re- 
places g with — g and adjusts the partitioning of pd so 



that atom A is anionic. Requiring continuity in the en- 
ergy at g = dictates that —SN*^ ^ — > 7(g) as g ^ 0^. 
Eq. (28) remains the same structurally. Equivalcntly, one 
could replace the subscript "A" with the subscript "B" 
everywhere in the procedure. The potential over the en- 
tire range of g is then represented in a piecewise fashion. 
While not as rigorous as a 3-state model, the present 
treatment does retain substantially greater simplicity. 

Another possible extension is to apply the above pro- 
cedure when the two resonance states both correspond 
to charged species. For CaO, for example, the effective 
g on Ca near equilibrium would be almost -1-2. As the 
CaO bond is stretched, g would decrease until it passed 
through a region of undetermined length where it would 
range between -|-1 and 0. Then there would be one form 
of Eq. (28) covering the range in g between and +1 
and a second form covering the range between +1 and 
+2 range. 

A final extension of this model would allow the entity 
B correspond to a more general environment than just 
one other atom. Most of the model presented here does 
not explicitly invoke the specific properties of a diatomic 
model. However, this extension requires separate consid- 
erations not pursued here. 



IV. DEFINITION OF PSEUDO-ATOM 
ENERGIES 

We now show how to define pseudo-atom energies for 
the 2-state model of the previous section. We do this by 
requiring consistency between the energies based on the 
density decompositions that were assumed in the previ- 
ous Section and the energies that would result from the 
corresponding wavefunction expressions. As with p, we 
assume that there is a decomposition F into pseudo-atom 
density matrices. In conformance with Rychlewski and 
Parr,^^ the decomposition applies to F rather than to H. 
Thus, for 

r = Fl + F^ , (29) 

the total energy decomposes into 

E = (H, Fl) + {H, F^) ^El + E^. (30) 

Analogous expressions are assumed to exist for each of 
the covalent and ionic contributions to the pseudo-atom 
energies. 

The pseudo-atom densities p\ and correspond to 
the pseudo-atom density matrices F^ and Fg. In the 
following relations, the expressions for atoms A and B 
arc analogous. Only the expressions for atom A will be 
given. We want the energies of each pseudo-atom defined 
via density matrices to be equal to the energies defined 
via densities. Consequently, we require that 

EI^EaIpI]. (31) 
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This identification places a new constraint on F^. Rig- 
orously speaking, it should be optimal in the sense of 
Levy CS-DFT .41.53 The constraint is that should yield 
the lowest energy for all other ensemble density matrices 
which integrate to p\: 



E' 



f ^ (i/, ri) = ^min, (i/, Ta) . (32) 



In terms of valence bond resonance states 



(33) 



Similar to Ep^, the pure-state terms in Eqs. (33), H*^ 
and -ffj* A, can be represented in conventional DFT lan- 
guage. The resonance energy cannot. Eq. (33) repre- 
sents pseudo-atom energies for any q. To evaluate the 
resonance energy, we again follow the EVB strategy of 
determining them from some particular value of q = go- 
The result is 



H. 



ci,A ~ 



e*pMq) - acc(go) g*c,A ~ «»»(go) HI 



2aci(go) 



si, A 



(34) 



If (7o is set equal to the optimum q for each value of atomic 
separation, then £'(go) = i?^(go) + £'b(9o) corre- 
spond to the experimental potential energy for AB. Note 
that all of the energies in traditional EVB are K val- 
ues. Some variations of EVB incorporate temperature- 
dependent solvent effects. 5^ We do not include these vari- 
ations here. Finite temperatures are not required to es- 
tablish the model. On the other hand, there is nothing 
here that precludes extending the analysis to a finite tem- 
perature ensemble. ^'i^'33'34. 59 

Substituting Eq. (34) into Eq. (33), an expression for 
the pseudo-atom energy is achieved. For any g, 

E*tS(i) = (ac»(9)/acj(go))i^A(9o) 

+ (acc(g) - acc(go) acj(g)/aci(go)) iilc,K 

+ (««(<?) - a«('?o) OLai((i) I OLaiHo)) Hl^ p^ . 

(35) 

Eq. (35) has the same structure as Eqs. (28), but refers 
to an individual atom. All of the quantities on the right- 
hand sides of these two equations can be deduced from 
experiment and/or decomposition calculations on reso- 
nance states and overlaps. 

Using the present formulations in simulations of larger 
systems naturally invokes consideration of chemical po- 
tential equalization. One can obtain a statement of chem- 
ical potential equalization from Eq. (35). The total en- 
ergy is given in Eq. (30). If E is minimized with respect 
to q, then small deviations from gopt will not change the 
total energy to first order: 







dEiq) 



dq g=gopt 

dEliq) 



dq 



dE^{-q) 



i3=9opt 



dq 



'?=9opt 



(36) 



dEliq) 



dq 



dE^{-q) 



9=9opt 



9='?opt 



-Kq) , (37) 



where fi is the chemical potential. Charge balance re- 
quires that the charges on A and B be exactly opposite. 
The minus sign in front of /U comes from the fact that q 
is related to the negative of the change in the number of 
electrons on A. 



V. DISCUSSION AND EXAMPLES 

Here we introduce approximations consistent with the 
dissociation limits of AB, present a general definition of 
the energy for a given charge, and discuss II2, HE, and 
LiH as examples. 



A. Neglect of Differential Overlap Model 

It is insightful to introduce a concept of neglect of dif- 
ferential overlap between resonance states, pd = 0, which 
we will refer to as NDOL. This leads to simple analytical 
expressions whose behavior can be examined in detail. 
NDOL is to be distinguished from zero differential over- 
lap (ZD0),3^ which refers to overlap between orbitals on 
different atomic centers. In the NDOL approximation, 
Pci = 0.^^ Under special conditions, ZDO implies NDOL. 

In the NDOL approximation, SN*^ ^ and Sd are zero. 
From Eq. (26), the dependence of 7 on g becomes 

7(g) = ±^q/il-q) . (38) 
From Eq. (9), the energy dependence on g becomes 



^NDOL(^) = ifec-2v/g(l-g)|i?c 
+ g(IPl-EA^) . 



(39) 



The — |ffcj| construct ensures that iJ^DOL corresponds to 
the ground state. Here we have used the fact that, by def- 
inition of the ionization potential IP and electron affinity 
EA, Ha — Hcc = IPa ^ EAg. These atom- in- molecule 
quantities include some electrostatic contributions. This 
quantity is also called a "bond hardness". 

Thus, Eq. (38) is consistent with previous energy ex- 
pressions obtained with Coulson charges^'^'^s^ss g^j^j with 
the 2-state model of Nalewajski^^. In Coulson, Murrell 
et. al,^'^ and McWeeny^* a fraction of ionic character is 
defined instead of a charge. That fraction is identical 
to 7^/(1 -f- 7^) which is g in the NDOL approximation. 
Clearly, at large R, Eq. (8) in combination with Eq. (26) 
approaches Eq. (39), Eq. (39) becomes linear in g, and 
Ha, Hcc, and approach asymptotic values (in R) 
analogous to MM's Eq. (2.40).^'' The linear behavior in 
the asymptotic regime is consistent with the conclusions 
of PPLB. Even more importantly, Eq. (39) is expressly 
non-analytical (i.e., it cannot be expanded in a Taylor 
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series) about q = and q == 1, whereas Eq. (8) in com- 
bination with Eq. (26) is analytical at both points away 
from the NDOL hmit ((5iV*^A 7^ 0)- Coulson's Fig. 5.7 
is a graph of jiq) vs. q, which clearly shows the non- 
analytical behavior, although he did not comment on it.'^^ 
The non-analytical behavior of the energy as a function 
of q is also seen in GC-DFT.^ Eq. (39) embodies the EVB 
representation of that behavior. The non-analytical be- 
havior results in a derivative discontinuity in the energy 
as a function of charge at integer values of the charge. 
Perhaps because orthogonality between resonance states 
is assumed, Nalewajski does not comment on the behav- 
ior of the derivative at integer charges.'''^ Cioslowski and 
Stefanov do see their version of the NDOL limit as con- 
nected with Likewise, PPLB,^ Perdew,^^ and 
Cioslowski and Stefanov^^ note that the derivative dis- 
continuity disappears once the systems in the GC rep- 
resentation begin to interact significantly. Eq. (26) cm- 
bodies that behavior as well. Examples of these charge 
dependences are illustrated below. 



Starting from Eq. (39) , we evaluate Hd at the optimum 
NDOL q, 

«-4(i-^). m 

where e = 2Hci/{Hii — Hcc)- The ground state corre- 
sponds to the negative root of 7(g) in Eq. (38), since 
as Hci 0, q must also go to zero, assuming that 
Hii — Hcc > 0. For convenience, we call the ground state 
value of Eq. (39) E^^^'^^ = i;NDOL(-^NDOL)^ The solution 
for the resonance energy is equivalent to the well-known 
EVB expression^^ 

\Hc^\ = ^(Hcc-EfJ^OL)^H,,-E^J^OL) . 

(41) 

Substituting Eq. (41) into Eq. (39), we find 



i?NDOL(^) = (1 - <Z) Hcc - 2 ^{1 - q){Hcc - E^POL) q(^Hu - E^J^OL) + q Hu . (42) 

I 



The dependence of E^^°^[q) on Hcc and Hu appears to 
be different from that implied by Eq. (28). In fact, by 
setting E{q^^'^^) = E^J^'^^ in Eq. (28), the two expres- 
sions become identical. 

We can gain further insight from Eq. (42) by solving 
for Hcc and H^.^ in terms of w Egs, i^xs, and 

Qopt ~ Qgs- In a typical diatomic, E^s corresponds to 
the first electronic state which dissociates to the ions A+ 
and B~. Because of the NDOL approximation, these 
three pieces of information are sufficient to specify the 
energy. To achieve the desired result, we first use the 
fact that Hci can be derived from either eigenvalue to 
find that Ha = E^s + Egg — Hcc- Next we solve for Hcc 
from Eq. (40) assuming that we know qg^. The result is 
that 

Hcc = {l- ggs ) Egg + qgg Ey,g (43) 

and 

Hu = qgs Egs + {1- qgg) E^^ . (44) 

In the NDOL approximation, the pure state energies are 
simple linear combinations of the cigcnenergies. Making 
all of the necessary substitutions and rearrangements in 
Eq. (42), we achieve the ensemble representation: 

i?NDOL(^) ^ ^ ^(^^ g^^) (^^^ „ ^ (45) 

where the occupation number is 

^{q, Qgs) = (Zgs - 2 ^q{l - q)qgs{l - qgs) + 9 (1 - 2ggs) . 

(46) 



To see that Eq. (45) has the desired properties, note that 
uj{q,qgs) lies between and 1 over the interval [0,1] in q, 
it is at g = qgs, and duj{q,qgs)/dq at q = qgs is also 
0. These properties are illustrated in Fig. 2. Taking into 
account our restriction to a 2-state model, if we apply 
atom decomposition to the eigenenergies, we would ob- 
tain the same form as Eq. (4), with uj^o = 1 — uj{q,qgs), 
lua+ = ^^iliQgs), and uj^- = 0. Importantly, the repre- 
sentation is in terms of eigenenergies instead of energy 
matrix elements. Furthermore, as qgs approaches 0, as 
in the assumed dissociation limit for AB, E^^'^^{q) be- 
comes linear in q. 

Another representation of the occupation number, 
Eq. (46), is significant. By using the relationship q = 
7^/(1 + 7^) from Eqs. (38), we obtain 

As is necessary physically, the occupation number is 
when 7 = 7gs, where 7gs corresponds to qgs- The com- 
plete generalization of Eq. (47) is equivalent to following 
the process steps outlined above. 

Eqs. (43) and (44) can be inverted. Inversion gives 
£;NDOL(^) in terms of Hcc and Hu: 

^NDOL ^ jj^^ ^ ^ (48) 

where 

Wcc = (1 - - 9gs)/(l - 2 qgs) (49) 
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FIG. 2: 2-state EVB-NDOL occupation number, Eq. (46), in 
the ensemble representation of E^^'~'^{q), Eq. (45). Charges 
are dimensionless. 



FIG. 3: 2-state EVB-NDOL coefficient for the covalent state 
given by Eq. (49). The ionic state coefficient is obtained by 
reflecting this Figure through cOcc ~ 1/2. Charges are dimen- 
sionless. 



and 



uju = {uj- ggs)/(l - 2gfgs) 



(50) 



Consistent with the conclusions of MM and physical 
necessity, the coefficients are not positive-semidefinite. 
Fig. 3 illustrates Eq. (49), the coefficient for the covalent 
state. In order to cover the range of energies between Egs 
and -Exs, the coefficients of Hcc and Ha cannot possibly 
be positive-semidefinite. The coefficients are not defined 
for Qgs = 1/2. At that value of the ground-state charge, 
Hcc niust equal Ha. 

It is interesting to see how the present results connect 
with classical electrostatic potentials (Eq. (3)). We can 
expand Eq. (46) locally as a function of q as long as the 
expansion point is not zero or one. First, 



dq 



The obvious value about which to expand is ^gg. By 
construction, duj{q,qgs)/dq\q^^ ~ 0. Second, 



d'^uj{q,qgs) 



I ggs (1 - ggs) 

2(g(l-g))3 



(52) 



If we evaluate Eq. (52) at q — ^gg, as (/gg ^ or 1, the 
expansion behaves badly. Alternatively, we might try 
expanding about q = 1/2, for the physically appealing 



reason that dE^^'^^{q)/dq\q=i/2 = IPl - EAg. Then 



d'^u{q,qgs) 



y'2gg3(l-ggg) . (53) 



9=1/2 



To second order, 

+ ((/gg - (1- a/2/8) y^ggg(l- ggs) 
+ (1 - 2 ggg - 9gg (1 - (7gg)/2) q 

+ ^2ggg(l-<7gg) <?V2) (iJxg - iJgg ) . 

(54) 

This quadratic expansion has the form of Eq. (3) and be- 
haves well physically under dissociation, even to the ex- 
tent of preserving the dissociation limit. The electroneg- 
ativity, hardness, and electrostatic contributions are em- 
bedded in the eigenenergies. The equivalent expansion of 
Eq. (48) might be more revealing in displaying these con- 
tributions. A key observation is that the quadratic term 
vanishes completely at infinite separation [q^s 0) and 
only the linear dependence survives. Again, the survival 
of the linear charge dependence is consistent with PPLB. 
There is no residual atomic hardness contribution as has 
appeared in many implementations of Eq. (3).^"^'^^'^* We 
speculate that a hardness contribution might be missing 
because we have considered only a 2-state model instead 
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of a 3-state model. Our reasoning behind this specula- 
tion is discussed below. We have ruled out the possibil- 
ity that introducing the NDOL approximation prior to 
making the expansion is a factor. Nevertheless, adding 
a third state must not change the fact that the coeffi- 
cient of the quadratic term must go zero in order for the 
results to be consistent with PPLB.^'^^''^^ Likewise, any 
hardness contributions to the linear term must be scaled 
by a coefficient that goes to zero at large R. Clearly, one 
cannot approximate either Eq. (26) or (39) quadratically 
to arbritrary accuracy in a global sense. 



is the reason behind our speculation that the isolated- 
atom hardnesses do not appear in Eq. (54) because of the 
limitations of the 2-state model. One advantage of the 
charge-generalized CS-DFT approach over a maximum 
entropy principle is that the true ground state proper- 
ties can be preserved in the same way that the ground 
state energy can be recovered — by minimizing over all 
allowed densities in conventional CS-DFT. 



C. Examples: HF, LiH, and Ha 



B. General Definition of the Energy 
for a Given Charge 

Next we consider the densities and their energies in the 
NDOL approximation. From Eq. (20) we have 



p''°°"(7) = (Pcc + 7'V.0/(l+7^) 



Thus, from Eq. (38) 

NDOL 



(?) = Pcc + q {Pii - Pec) 



(55) 



(56) 



By symmetry in Eq. (55), ^(7) = p(— 7), but, from 
Eq. (8), E{-i) ^ E{~-i). To address this difficulty, 
MM appeal to a maximum entropy principle. We prefer 
instead to appeal to CS-DFT.^^'^i.ss Accordingly, CS- 
DFT instructs us to place the r'^°'^^(7) into groups de- 
fined by the density that they produce. Since r'^'^'^'^(7) 
and rNDOL(_^) i^Q^i^ yjgi^j ^jjg gg^j^g density, p^^'^^iq), 

they are grouped together. The energy assigned to 
^NDOLj^NDOL^^-jj is the minimum of the energies for 
these two F's: 



E 



NDOLr^NDOL 



NDOL 



(7(<Z))), 



(ij,r 



NDOL 



(-7('?))>} . (57) 



The ground state energy is then the minimum over all q 
of ^NDOL[pNDOL(^)]^ When differential overlap (p„ ^ 0) 
is included, the densities appear to become unique for 7's 
of different signs. 

Note that a similar procedure could be followed for 
MM's 3-state case at the ZDO level that they assume. 
See also Ref. 33. Even with this comparatively simple 
extension, the situation is less clear than in the 2-state 
case. Many densities may have the same charge on the 
atoms. By again appealing to CS-DFT, one can assign 
E{q) for a given q by minimizing over the energies of all 
densities with the same q. That is, as a straightforward 
extension of CS-DFT, we very generally define 



E{q) 



mm 



E 



CS-DFT 



[p{q,S)] 



where 5 represents all of the other undetermined parame- 
ters of the density p{q, 6). By minimizing over <5, one may 
introduce dependencies on the energy matrix elements 
that arc absent from the present 2-state model. ^'^ This 



First we discuss modeling the polar molecules HF and 
LiH in the NDOL approximation. Then, for the non- 
polar molecule H2, wc compare the NDOL and general 
cases. There we utilize the Weinbaum wavefunction.^^'^^ 
It provides an excellent illustration of the ambiguities en- 
coountered in defining the resonance state wavefunctions. 
We also examine various approximations for 5N*^ ^ and 
examine the change in the charge dependence as a func- 
tion of R. 

For HF and LiH, wc use RKR curves^^ to define E^^ 
and E'xs. For the X^Y.+ and B^Y.+ states of HF, the RKR 
data are from Di Lonardo and Douglas, for the X^I]+ 
and A^I]+ states of LiH, the RKR data are from Chan et 
alP and Pardo et alJ^ We use the calculations of Ref. 19 
to define qgs- These data are shown in Figs. 4 and 5. The 
B^I]+ and A^E+ states dissociate to ions, H+ and F~ 
and Li"*" and H~ , respectively. To allow matching of the 
different spatial ranges of the data, analytical fits for the 
energy curves were made with the Rose and Rose-|-ionic 
functional forms. The charge data were fit with the 
functional form, Ko + {Ki- Ko){q^'' / {1 + {q^^ + iff" ) ) ) , 
where the K^s are fitting parameters. Hcc and Ha were 
computed from Eqs. (43) and (44), respectively. In both 
cases, Hcc and Ha cross at qgs = 1/2 and meet their 
respective states at the dissociation limits. Representa- 
tive shapes of the charge dependence for given values of 
R are shown in Figs. 6 and 7. The correct dissociation- 
limit behavior is observed in both cases. In simulations 
where each atom remains within a unit charge interval, 
the NDOL model might therefore prove to be useful, al- 
though it is unlikely to be quantitative. The difference 
in charge transfer characteristics between the present re- 
sults and PPLB can be seen by contrasting the i? = 3 A 
curve of Fig. 7 with the i?c = 3.1 A curve of their Fig. 1.^ 
The transfer is more gradual here and passes through 
fractional charge states, compared to PPLB which is 
very sharp and passes directly from a completely covalent 
state to a completely ionic one. 

Next we consider H2. The simplest valence bond 
(58) form for the covalent state ip'c'^ ~ (0a(1)'/'b(2) -I- 
(/)b(1)0a(2))/(2 + 2S'iB)l/^ where S'ab is the atomic 
orbital overlap.^" Our naive inclination for the ionic 
state is to use the familiar = ((/)a(1)0a(2) + 

0b(1)0b(2))/(2 + 2S'iB)i/2^ The total wavefunction is 
ip = cijjjc + l^i)- ft turns out that the total densities for 



11 




R(A) 



6 - 



> 

a 
c 

UJ 



-4 - 



R(A) 

- 1 

- 2 

- 3 

. 4 



FIG. 4: EVB model of HF. All energies are in eV, distances 
in A, and charges dimensionless. RKR data from Ref. 76. 
Charges from Ref. 19. 



0.0 




R(A) 



FIG. 5: EVB model of HF. All energies are in eV, distances 
in A, and charges dimensionless. RKR data from Ref. 77 and 
78. Charges from Ref. 19. 
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FIG. 6: Eq. (45) at discrete values of R for HF. All energies 
are in eV, distances in A, and charges dimensionless. 
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both states are identical: pcc = Pa = (0a + 2 Sab 4>a<Pb + 
'^b)/(1 + '^ab)- (Designation of the electronic coordinate 
is suppressed in the densities for readability.) The inter- 
ference density is pd = ((0a + '?^b) '^'ab + 2 0a0b)/(1 + 
Sab)- 

If we assume a simple Hirshfeld partitioning as our to- 
tal density decomposition strategy/^'^^'^^ then it is natu- 
ral to assume that each component of the density is like- 
wise scaled by 4>\I{'P\ + 06) order to obtain the atom 
A contribution. However, we find that because of symme- 
try, 5N*^ A = 0; which cannot be correct. Of course, the 



FIG. 7: Eq. (45) at discrete values of R for LiH. All energies 
are in eV, distances in A, and charges dimensionless. 



problem is that the assumption about the Hirshfeld par- 
titioning form is incorrect. This partitioning cannot lead 
to a value of pu^A that integrates to Nao — 1 = 0. A dif- 
ferent partitioning must be chosen to force all of the den- 
sity in Pa to belong to atom B. In addition, we know the 
optimum values of 7 for these two resonance states^^'^^ 
and that these values correspond to q = Q. For instance. 
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at the equilibrium H2 separation, 7opt = —0.26, with the 
resonance energy defined as positive. From Eq. (23), 
we also can deduce that SN*^ ^ = 7opt/2 = —0.13. This 
analysis hints at the subtle properties that the density 
decomposition must possess if one implements explicit 
decomposition of pd- 

Another immediate insight is that any fixed value of 
7 applied within ijj leads to an acceptable covalent state. 
Clearly, the energy change caused by a polarization of the 
H2 charge density will differ depending on one's choice for 
the resonance state. Similarly, the density decomposition 
will show some sensitivity to this choice. 

Of course, a more sensible choice for the ionic state 

(2.) 

is ipl = 0b(1)0b(2) because, for a homonuclear di- 
atomic, the charged states correspond to broken charge- 
symmetry states. Even in this simple case, the partition- 
ing of pci may be nontrivial. 

As a final consideration, we illustrate the influence of 
overlap on ionicity. For this purpose, we utilize our orig- 
inal choices for covalent and ionic wavefunctions, t/'c^'' 
and ipi^^ ■ The overlap integral is estimated in terms of 
Sab as given above. Taking 0(r) = \/a^/7rexp(— ar), 
the atomic overlap is S'AB(a) = exp(— a)(l -I- a -I- 1/3 a^), 
where a = aR. The effect of overlap is shown in Fig. 8 
as 7('i')(l — (?). There it can been seen that the NDOL 
approximation becomes accurate beyond approximately 
2 A. The relatively large range of R for which the NDOL 
approximation is accurate in this case is due in no small 
part to the fact that Sd depends on the square of Sab- 
This broad range of accuracy and the fact that Eqs. (45) 
and (46) are well-behaved for any R makes it tempting to 
use NDOL in general. However, near equilibrium bond 
lengths. Fig. 8 indicates that overlap effects should not 
be ignored. 



VI. CONCLUSION 

In order to apply chemical potential equalization in 
a simulation that involves conditions far from reference 
states, the potential energy must be defined for arbi- 
trary values of the charges. To address this challenge, 
we have derived a new charge-dependent pair potential 
from a 2-state empirical valence bond model. The charge 
is defined from a decomposition of the density into con- 
stituent contributions. The explicit charge dependence 
is deduced by requiring consistency between the density 
decomposition and the wavefunction descriptions of the 
ground state. The energy expression can be made valid 
for any range of charge of interest. The decomposition 
theme is further extended to define the energy of indi- 
vidual constituents, again by requiring consistency be- 
tween the density decomposition and the wavefunction 
descriptions. The energy of the system for a given value 
of charge is made unique by appealing to constrained 
search density functional theory. An examination of the 
model shows linear dependence on charge at the disso- 




q (dimensionless) 

FIG. 8: Ionicity scaled by 1 — g at discrete values of R for H2. 
Distances are in A and charges are dimensionless. 



ciation limit, as well as discontinuous behavior in the 
derivative of the energy as a function of charge at inte- 
ger values of the charge. This behavior is consistent with 
the analysis of Perdew et al.^ and others. ^^■'^^ The po- 
tential energy for arbitrary charge and separation of the 
constituents can be represented as an ensemble average of 
the eigenenergies with a nonlinear, analytical dependence 
of the occupation number on charge. The representation 
of the potential energy in terms of the pure-state energy 
matrix elements is possible, but, by physical necessity, 
cannot be expressed as an ensemble average with posi- 
tive semidefinite coefficients. 

To determine a pair potential for all values of q and R 
with this method, one needs five reference or calibration 
curves. These most often will be computed values of the 
ground-state energy, the charges along the ground-state 
energy curve, the resonance state overlap integral, the co- 
valent energy and ionic energy. Measurements of ground 
and ionic-excited state energies and charges can also be 
used. The covalent energy and ionic energy are deter- 
mined for integer charges only. To determine a point on 
the potential energy curve, a well-defined five step proce- 
dure is followed using the five input curves. An immedi- 
ate application of the method could be to construct the 
reference potential curves for the A2, B2, and AB sys- 
tems. Simulations on arbitrary mixtures of these three 
types of systems under nonequilibrium initial conditions 
would then be possible. Chemical potential equalization 
would be used to dynamically adjust the charges of the 
constituents. More generally, we envision the present ap- 
proach as forming the basis for a new class of charge- 
dependent empirical potentials for use in large-scale sim- 
ulations of reactive systems. 



13 



Acknowledgments 

The work of S.M.V. was performed in part at Los 
Alamos National Laboratory under the auspices of the 
U. S. Department of Energy, under contract No. W- 
7405-ENG-36, and funded through its Center for Semi- 
conductor Modeling and Simulation, a CRADA program 
performed jointly with the Semiconductor Research Cor- 



poration, and through the Advanced Fuel Cycle Initia- 
tive. S.M.V. thanks the University of New Mexico, De- 
partment of Physics and Astronomy for its hospitality 
during the 2003-2004 academic year. S.R.A. would like 
to thank the National Science Foundation for support 
during the initial stages of this work (DMR-9520371). 
This work was supported by National Science Founda- 
tion grant No. CHE-0304710. 



^ R. T. Sanderson, Science 114, 670 (1951). 

^ J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, Jr., 

Phys. Rev. Lett. 49, 1691 (1982), hereafter referred to as 

PPLB. 

^ R. G. Parr and R. G. Pearson, J. Am. Chem. Soc. 105, 
7512 (1983). 

" S. W. Rick, S. J. Stuart, and B. J. Berne, 

J. Chem. Phys. 101, 6141 (1994). 
^ A. K. Rappe and W. A. Goddard, III, J. Chem. Phys. 101, 

6141 (1994). 

F. H. Streitz and J. W. Mintmire, Thin Solid Films 253, 
179 (1994); Phys. Rev. B 50, 11996 (1994); Langmuir 12, 
4605 (1996). 

D. T. Nguyen, A. C. Scheiner, J. W. Andzelm, 
S. Sirois, D. R. Salahub, and A. T. Hagler, J. Com- 
put.Chem. 18, 1609 (1997); L. F. Pacios and P. C. Gomez, 
J. Molec. Struc. (Theochem) 544, 237 (2001). 

* S.-Y. Sheu, D.-Y. Yang, H. L. Selzle, and E. W. Schlag, 
J. Phys. Chem. A 106, 9390 (2002). 

^ S. Hammes-Schiffer, Acc. Chem. Res. 34, 273 (2001). 
M. L Page, "The mechanism of chemical catalysis used 
by enzymes," in New Comprehensive Biochemistry, Vol. 6 
The Chemistry of Enzyme Action edited by M. I. Page, 
(Elsevier, Amsterdam, 1984), pp. 229-270. 
S. J. Benkovic and S. Hammes-Schiffer, Science 301, 1196 
(2003). 

J.-W. van der Horst, P. A. Bobbert, P. H. L. de 
Jong, M. A. J. Michels, G. Brock, and P. J. Kelly, 
Phys. Rev. B 61, 15817 (2000); M. Rohlfing and 
S. G. Louie, Phys. Rev. B 62, 4927 (2000); M. Rohlfing, 
M. L. Tiago, and S. G. Louie, Synth. Metals 116, 101 
(2001). 

R. G. Parr and L. J. Bartolotti, J. Am. Chem. Soc. 104, 
3801 (1982). 

" R. F. Nalewajski, J. Am. Chem. Soc. 106, 944 (1984). 
R. G. Pearson, Hard and Soft Acids and Bases, (Dowden, 
Hutchinson, and Ross, Stroudenberg, PA, 1973). 
R. F. Nalewajski and R. G. Parr, J. Chem. Phys. 77, 399 
(1982). 

R. F. Nalewajski and M. Koniriski, J. Phys. Chem. 88, 
6234 (1984). 

F. L. Hirshfeld, Theor. Chim. Acta 44, 129 (1977). 

J. Cioslowski and B. B. Stefanov, J. Chem. Phys. 99, 5151 

(1993). 

^" U. W. Schmitt, and G. A. Voth, J. Phys. Chem. B 
102, 5547 (1998); U. W. Schmitt, and G. A. Voth, 
J. Chem. Phys. Ill, 9361 (1999); M. Cuma, 
U. W. Schmitt, and G. A. Voth, J. Phys. Chem. A 
105, 2814 (2001). 

A. Alavi, L. J. Alvarez, S. R. Elliott, and L R. McDonald, 



Phil. Mag. B 65, 489 (1992). 
^■^ B. W. H. van Beest, G. J. Kramer, and R. A. van Santen, 
Phys. Rev. Lett. 64, 1955 (1990). 

Y.-P. Liu, K. Kim, B. J. Berne, R. A. Friesner, and 
S. W. Rick, J. Chem. Phys. 108, 4739 (1998). 
R. P. Iczkowski and J. L. Margrave, J. Am. Chem. Soc. 83, 
3547 (1961). 

G. Klopman, J. Chem. Phys. 43, S124 (1965). 

R. G. Parr and W. Yang, Density Functional Theory of 

Atoms and Molecules (Oxford, New York, 1989). 

L. von Szentpaly, J. Mol. Struct. (THEOCHEM) 233, 71 

(1991). 

L. von Szentpaly, Chem. Phys. Lett. 245, 209 (1995). 
^® L. von Szentpaly and D. O. Niel Gardner, J. Phys. Chem. A 

105, 9467 (2001). 
^° W. J. Mortier, S. K. Ghosh, and S. Shankar, 

J. Am. Chem. Soc. 108, 4315 (1986). 

P. Geerlings, F. De Proft, and W. Langennaeker, 
Chem. Rev. 103, 1793 (2003). 

L. Jaroszewski, B. Lesyng, and J. A. McCammon, 
J. Mol. Struct. (THEOCHEM) 283, 57 (1993); P. Gro- 
chowski, B. Lesyng, P. Bala, and J. A. McCammon, 
Int. J. Quant. Chem. 60, 1143 (1996); J. Trylska, P. Gro- 
chowski, and M. Geller, Int. J. Quant. Chem. 82, 86 
(2001). 

R. F. Nalewajski, Int. J. Quant. Chem. 69, 591 (1998). 
J. Morales and T. J. Martinez, J. Phys. Chem. A 105, 
2842 (2001). 

J. P. Perdew, in Density Functional Methods in Physics, 
NATO Advanced Science Institute Series, Vol. 123, edited 
by R. M. Dreizler and J. da Providencia (Plenum Press, 
New York, 1984). 

C. A. Coulson, Valence, 2nd edition (Oxford University 
Press, Oxford, 1961). 

J. N. Murrell, S. F. A. Kettle, and J. M. Tedder, Valence 
Theory (John Wiley and Sons, 1965) . 

R. McWeeny, Coulson's Valence (Oxford University Press, 
Oxford, 1979). 

P. Hohenberg and W. Kohn, Phys. Rev. B 136, 864 (1964). 
^° W. Kohn and L. Sham, Phys. Rev. 140, A1133 (1965). 
'^^ M. Levy, Proc. Natl. Acad. Sci. USA 76, 6062 (1979). 
'^^ W. Mofiit, Proc. Roy. Soc. (London) A210, 245 (1951). 
■''^ R. F. W. Bader, Atoms in Molecules: A Quantum Theory 

(Oxford University Press, Oxford, 1990). 
'''' R. G. Parr, R. A. Donnelly, M. Levy, and W. E. Palke, 

J. Chem. Phys. 88, 3801 (1978). 

M. F. Guse, J. Chem. Phys. 75, 828 (1981). 

L. Li and R. G. Parr, J. Chem. Phys. 84, 1704 (1986). 

M. J. Mehl, L. L. Boyer, and H. T. Stokes, J. Phys. 

Chem. Solids 57, 1405 (1996); H. T. Stokes, L. L. Boyer, 

and M. J. Mehl, Phys. Rev. B 54, 7729 (1996). 



14 



A preliminary version of this work has been presented. 

S. R. Atlas and S. M. Valone, Bull. Am. Phys. Soc. 47, 

Part II, 1213 (2002); S. R. Atlas and S. M. Valone, to be 

submitted (2003). 
'^^ P.-O. Lowdin, J. Chem. Phys. 18, 365 (1950). 
^° R. S. MuUiken, J. Chem. Phys. 23, 1833 (1955). 
5^ C. A. Coulson and U. Danielsson, Ark. Fys. 8, 239 (1954). 

A. Warshel and A. Bromberg, J. Chem. Phys. 52, 
1262 (1970); A. Warshel and R. M. Wiess, J. Am. 
Chem. Soc. 102, 6218 (1980); J. Aqvist and A. Warshel, 
Chem. Rev. 93, 2523 (1993). 

^'■^ S. M. Valone, J. Chem. Phys. 73, 1344 (1980); 
J. Chem. Phys. 73, 4653 (1980). 

R. S. MuUiken, Phys. Rev. 50, 1017 (1936); ibid., 50, 1028 
(1936). 

J. Q. Broughton and M. J. Mehl, Phys. Rev. B 59, 9260 
(1999). 

S. Weinbaum, J. Chem. Phys. 1, 593 (1933). 

C. A. Coulson and I. H. Fischer, Phil. Mag. 40, 386 (1949). 

B. Barbiellini and A. Shukla, Phys. Rev. B 66, 235101 
(2002); T. K. Chanty, V. N. Staroverov, P. R. Koern, and 
E. R. Davidson, J. Am. Chem. Soc. 122, 1210 (2000). 

A. Cedillo, P. K. Chattaraj, and R. G. Parr, 

Int. J. Quant. Chem. 77, 403 (2000). 
®° K. Ruedenberg, Rev. Mod. Phys. 34, 326 (1962). 

R. F. Nalewajski and R. G. Parr, Proc. Natl. Acad. 

Sci. USA 97, 8879 (2000). 
®^ R. F. Nalewajski and R. Loska, Theor. Chem. Acc. 105, 

374 (2001). 

^'^ R. F. Nalewajski and R. C. Parr, J. Phys. Chem. A 105, 
7391 (2001). 

J. W. Storer, D. J. Giesen, C. J. Cramer, and D. G. Truh- 
lar, J. Comput.-Aided Molec. Design 9, 87 (1995). 
"■"^ J. Li, T. Zhu, C. J. Cramer, and D. G. Truhlar, 
J. Phys. Chem. A 102, 1820 (1998). 



Q. Zhao, R. C. Morrison, and R. G. Parr, Phys. Rev. A 

50, 2138 (1994). 
•^■^ P. W. Ayers, J. Chem. Phys. 113, 10886 (2000). 
®* Y. Wang and R. G. Parr, Phys. Rev. A 47, R1591 (1993). 
®^ M. I. Baskes, Phys. Rev. Lett. 59, 2666 (1987). 
™ M. I. Baskes, J. S. Nelson, and A. F. Wright, Phys. Rev. B 

40, 6085 (1989). 

M. I. Baskes, Phys. Rev. B 46, 2727 (1992). 

J. Rychlewski and R. G. Parr, J. Chem. Phys. 84, 1696 

(1986). 

J. M. Parks and R. G. Parr, J. Chem. Phys. 28, 335 (1958). 
For an arbitrary set of basis wavefunctions, there is no sim- 
ple way to determine which wavefunctions map to which 
densities. Even in a finite basis set, pathological cases are 
possible. 

"^^ R. Rydberg, Z. Physik 73, 376 (1931); O. Klein, Z. Physik 
76, 266 (1932); A. Rees, Proc. Phys. Soc. (London) A 59, 
998 (1947). 

''^ G. Di Lonardo and A. E. Douglas, Can. J. Phys. 51, 434 
(1973). 

■^■^ Y. C. Chan, D. R. Harding, and W. C. Stwalley, 
J. Chem. Phys. 85, 2436 (1986). 

A. Pardo, J. J. Camacho, and J. M. L. Poyato, 
Chem. Phys. 108, 15 (1986). 

J. H. Rose, J. Ferrante, and J. R. Smith, 
Phys. Rev. Lett. 47, 675 (1981); J. H. Rose, J. Ferrante, 
and J. R. Smith, Phys. Rev. B 28, 1835 (1983); J. H. Rose, 
J. R. Smith, F. Guinea, and J. Ferrante, Phys. Rev. B 29, 
2963 (1984); J. Ferrante, and J. R. Smith, Phys. Rev. B 
31, 3427 (1985); J. R. Smith, H. Schlosser, W. Leaf, 
J. Ferrante, and J. H. Rose, Phys. Rev. A 39, 514 (1989). 
*° J. C. Slater, Quantum Theory of Molecules and Solids, 
Vol. 1: International Series in Pure and Applied Physics, 
L. I. Schiff editor, (Mcgraw-HiU, New York, 1963). 



